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1. Introduction 

Presently, there is much interest in the problem of stability of matter confined by its self- generated gravity. This 
self -consistent approach dates back to the geons of Wheeler [1]. Recently, the work of Lee et al. [2,3] stimulated 
further progress. They pointed out that a star, regarded as a gravitational soliton [4 — 6], can have a mass which is 
larger than the Chandrasekhar type limit for gravitational collapse [7, 8, 9]. This opens up a new avenue for studying 
the structure of a star under unusual matter conditions. Its stability is the most important question. So far, the 
dynamical stability of boson stars has been analysed [10 — 12, 13] by means of perturbation theory. In this paper, 
we will apply a method which was proposed by one of us ([14] and references therein) for nongravitational solitons. 
In general, the method consists of investigating the critical points of a mapping and the construction of bifurcation 
diagrams. In our case, a two-dimensional subspace of the dynamical variables of the boson field is mapped into the 
space of the integrals of motion, such as the gravitational mass M and the total particle number N. Using Arnold's 
classification [15] of singularities of differential maps (catastrophe theory), we are able to derive general criteria for 
the stability of the star. 

2. Coupled Einstein— scalar field equation 

As a general-relativistic model of a boson star, we consider a self-interacting scalar field $ describing a state with 
zero temperature. This field is self-consistently coupled to its own gravitational field via the Lagrangian 

C = l-^\R + ^^\ [g^^id,'f*){d.<f)-Ui\-f\')] , (1) 

where n = SttG is the gravitational constant in natural units, g the determinant of the metric g^^, p.,v = (0, 1, 2, 3), 
R the curvature scalar, and U{\<^ ]^) the self-interaction potential. We will investigate to what extend the form of U 
infiuences the stability of the star. 

From the principle of least action we obtain the coupled Einstein-Klein-Gordon equations: 

R^.u - \gt,.R = -kT^u{^), (2) 
□+^)^ =0, (3) 

where 

T^j.vi'^) = {d^^*){dij<^) - {g,iu/ \/\ g |)'C($) is the energy-momentum tensor and □ = 
the generally covariant d'Alembertian. 

In this paper, we restrict ourselves to the static, spherical symmetric metric 

rfs2 = e^^'-Ut^ - e^^^'Ur'^ - r'^{de^ + sin^ Odcj)'^), (4) 

in which the functions v = z/(r) and A = A(r) depend on the Schwarzschild type radial coordinate r. For the boson 
field, we make the stationarity ansatz 

$(r, t) = P(r)e-''"* , (5) 

which describes a spherically symmetric bound state with frequency to. The resulting coupled system reads 

v' + X' =K,{p + Pr)re^ , (6) 
A' =Kpre^-ie^ + i , (7) 

P"+ Q(j.'_A') + ^^ P' = e^^P-e^-''u;2p . (8) 
The energy-momentum tensor becomes diagonal, i.e. T^^ = diag {p, —Pr, —p±, —p±) with 

p = i(w2p2g-. ^ p/2g-A ^ ^ (9) 
Pr =p-U ,p^=Pr- P'^e-^ . (10) 

The form of is familiar from an ideal fluid, except that the radial and tangential pressure generated by the scalar 
field are in general different, i.e. Pr=/=p±. This fractional anisotropy aj := (p^ ^ P±)/Pr has already been noted by 
Ruffini and Bonazzola [5]. Moreover, Gleiser [13] found that all boson stars have the same amount of anisotropy at 
the radius of the star. 
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Because of the contracted Bianchi identity V — ^g^j/i?) = 0, a further equation involving = tJ' = —p± 
is identically satisfied. Eq. (2.7) possesses a Schwarzschild-type solution 

r 

e-Hr) = 1 _ air) := J px^dx. (11) 



where a{r) is the mass function. For the polynomial self-interaction 

f/:=m2|$|2 + ia|$|^+i/3|$|^ (12) 

these equations has been solved numerically for nonsingular, finite mass and zero-node solution [3,4,5,16,17]. Two 
and higher node solutions occurred already in Rcf. [6]. For a massless scalar field with U = 0, an exact solution is 
known, cf. Ref. [18]. For a massless real scalar field, Christodoulou [19] could show that a spherically symmetric 
time-dependent field configuration must either disperse to infinity or, for non-vanishing Bondi mass, forms a black 
hole. 



3. Integrals of motion 

The concept of an energy-momentum 4-vector for a field configuration is a notoriously subtle [20] in general 
relativity. However, the exponential decrease of the radial function P(r) ~ exp[—^/m'^ — uP- r] for | w |< m yields an 
isolated, static system for which the Tolman mass formula 

oo 

M:= J (2To° - T/) vTffT = An J [2oj^P^e-'' - U] e^'+^^/^^dr (13) 



applies (cf. [6]). It can be derived from the local conservation law dviT^/ + t^^) = 0, where T^" = ^ and 
T^^ is the gravitational energy-momentum complex. For a boson star, the explicit expression (13) does not involve 

derivatives, in contrast to the Schwarzschild mass MschwarzschUd '■= 47rQ:(oo) that is commonly studied [11] in this 
context. Friedberg et al. implicitly rederived the equivalence of the Tolman and the Schwarzschild mass (see (2.27) 
of Ref. [3a]; cf. also [7]). 

A second "integral of motion" arises from the fact that the Lagrangian (1) is invariant under the global phase 
transformation <l> $e~*''. Therefore the Noether current density 

f = ^WaJg'"' [^*d.^ - $a,$*] (i4) 

is locally conserved, i.e. = 0. The time-component j° integrated over space yields the particle number N or the 

charge Q: 

•DC' 

N=9.=4wco f e<^^-''^/\^P^dr . (15) 



Since the current density (14) is a "measure" for the radial distribution of the "particles" in the boson star, its 
effective radius can be defined by 



i? := 1 1 ri'^dS, = e('-'^y'r'P'dr 



(16) 



On account of the fractional anisotropy a/ another interesting radius Rq could be obtained from the node p±{Ro) = 
in the tangential pressure p± . This radius Rq separates the interior part of the boson star from a marginal layer in 
which p± becomes negative before it decreases exponentially [17]. 
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4. Smooth mapping (Whitney surface) 

In order to investigate the stability of soUton-type solutions against radial perturbations, we consider the two- 
dimensional mapping 

F : (k,uj) ^ {M,N) , (17) 

where A: is a variational parameter which dilatates the radius R of the star and lo the frequency eigenvalue. The 
parameter k induces a scaling of the metric, the frequency, and the scalar field in accordance with their normal 
physical dimensions 

ds"^ k^ds^ , uj w/k, P{r) kP{kr) , (18) 

such that the particle number N is kept fixed. 

In order to classify the singularities of this mapping F, let us consider the Jacobi matrix 

fdM/dk dM/duj\ , . 

\ dN/dk dN/du) ) ' ^ > 

According to Whitney's theorem [14. 15], the singularities of the mapping F can be one of three types, depending on 
the rank Rj = 2, 1 and 0, respectively. Since we require the soliton solution to be an extremal point of the Lagrange 
manifold, we have: 

dM ON 

-=0 - = 0. (4.2) 

For the soliton the rank of J is Rj < 2 and, consequently, the singularities of the mapping F may have cither Rj = 1 
or Rj = 0. In that case, our soliton solution corresponds to the extremal or critical points of the Whitney surfaces 
which has a very definite form (see [14, 15]). In numerical examples, the dependence 

M = M{u) , N = N{uj) , with w = w(cr(0)), (4.3) 

on the frequency lo can be smoothly converted into a function of the central density (t(0) = y^P(O) such that the 
critical points coincide (Fig. 1). If the rank of J is zero, the critical points arc degenerate. The maxima and minima 
of M = M(ct(0)) and = N{a{Q)), see Fig. 1, correspond to the A2 singularity, in the notation of Arnold. Other 
points of the curves in Fig. 1 correspond to the critical points Ax . 



5. Bifurcation diagram 

In order to classify the nondegenerate Ai, we need to consider the bifurcation diagram M = M{N) [3]. Equivalently, 
we may consider the binding energy 

B = M-mN = B{N) (5.1) 

as a function of the particle number (Fig. 2). A further "magnification" of B{N) is achieved in Fig. 3. According 
to the Whitney theorem, the cuspoidal points of these diagrams classify the A2 singularity, whereas the other points 
of the diagram correspond to the Ai singularity. Each cusp represents some Whitney surface, which is a part of 
the mass-energy surface. As shown in Ref. [14], the minimum on this surface corresponds to the stable soliton, the 
maximum corresponds to the unstable soliton. At the cuspoidal point, the minimum coalesces with the maximum, and 
the soliton loses its stability. Thus, the lower branch of the lowest cusp corresponds to the absolutely stable soliton. 
The upper branch of the first cusp, which is, at the same time, the lower branch of the second cusp, corresponds 
to the unstable soliton. The upper branch of the second cusp also corresponds to unstable solitons, which, however, 
suffer from a difiierent kind of instability than the soliton of the lower branch of the second cusp. 

For the boson star, the degrees of freedom of the configuration space are very large. The fact that the lower branch 
of the lower cusp corresponds to absolutely stable soliton means that there we have minima for all directions in the 
configuration space. 

The higher branch of the first lower cusp corresponds to a maximum. This maximum occurs in that section of 
the mass-energy surface which depends on the radius R of the star. For the second cusp the appearance of a new 
instability depends on the mutual branching to other cusps. There are the two possibilities that the next (third) 
branch goes higher or lower than the second one. In the first case according to Whitney's theorem - the minimum 
transforms into a maximum at the transition from the second branch to the third one, and new instability appears. 
Vice versa, in the second case, one maximum transforms into minimum after the transition through the cuspoidal point 
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and one instability disappears. The numerical data (Fig. 3) show that for the third cuspoidal point one instability 
disappears, in accordance with the general picture developed here. This disappcarcncc of one instability has not been 
pointed out in previous works on the stability of boson stars. 

In the literature [3, 21] the corresponding diagrams are obtained by applying the method of small perturbations 
to the problem of stability. In our approach, the bifurcation diagram is a key point of the analysis and is gained 
by analyzing the topology of the Whitney surface. In the first step of the analysis, we calculate the bifurcation 
diagram for the mapping which relates the integrals of motion M and N to some degrees of freedom of the boson 
star. Following [14] we connect this bifurcation diagram to the mass-energy surface M = M{R, N). 

The mass-energy surface described by the bifurcation diagram corresponds to some complicated manifold of catas- 
trophe, which probably has never been seen before in the theory of singularities. The order of this grand catastrophe 
depends on the number of cusps in the bifurcation diagram. The static description of the bifurcation diagram is that 
at every cusp there is a transition from minimum to maximum or vice versa. For every cusp there occurs one Whitney 
surface. As an application to our case, we find at the first cusp that a minimum transforms to a maximum. This 
maximum will not be affected at the following cusps, if at the second cusp the following branch goes to higher mass 
values. If this branch reached lower mass values, we will get a more complicated "Whitney surface" in analogy to the 
case of the swallow's tail [22] . 

For neutron stars or white dwarfs, we get a complicated bifurcation diagram consisting of the many cusps (see Fig. 
8 of Ref. [8]). Again, each cusp corresponds to a Whitney surface. The first four lower branches of this bifurcation 
diagram (in the direction of increasing density) describe a section of manifold of catastrophe which is usually called 
a butterfly. Consequently, for these four branches this two-dimensional manifold of catastrophe corresponds to the 
two-dimensional Whitney surface in the bosonic case. 

For the second cusp of the bosonic bifurcation diagram there is a Whitney surface also, i.e. at a cuspoidal point of 
this cusp there is a transition from minimum to maximum. At the following cusps, if there is a transition to higher 
mass values, the maxima stay maxima. We understand such applications of catastrophe theory not only for the boson 
star, but also for fermion Q-stars, white dwarfs, and neutron stars. The latter have been considered by Harrison, et 
al. [8] . They 'saw' the Whitney surface without drawing knowledge from catastrophe theory. 

The simplest way to identify the instabilities, which have been qualitatively predicted by catastrohpe theory, is to 
consider perturbation theory. Because the perturbative equations set up a Sturm-Liouville eigenvalue problem, the 
characteristic frequencies of the perturbation series have increasing absolute values. This holds also in the case of the 
boson star [11, 12, 13]. Moreover, it was shown for the boson star [11] and for the neutron star [8] that, at each cusp, 
one of these frequencies changes sign. Thus, each such instability can be identified with a corresponding one obtained 
from applying catastrophe theory. 

Until now we have considered only the mapping of the two-dimensional space {k, w) into two-dimensional space 
{M,N). Since the first space counts the number of the degrees of freedom of the star, we may extend the domain 
(fc,w) of the mapping by including the characteristic frequencies. In the next step, new states of the star with 
different dynamical behavior could be taken into account than those we see already in the bifurcation diagram. 
These states correspond to other points on the mass-energy surface. For example, for burning stars oscillations with 
large amplitudes occur at their finite stage as red giants. The evolution of these oscillations cannot be described by 
perturbation theory. This kind of dynamical behavior of the star will bo discussed in Section 6. 

In fiat spacetime, the dependence M on N has been investigated by Friedberg et al. [21]. Although they pointed 
out that the minimal energy branch of M versus N is stable, which is true, they have not investigated the stability 
of these solitons for all values M and N. 

6. The different regimes of the stctr's behavior 

In our method, like in the theory of singularities of smooth mappings, the bifurcation diagram plays an important 
role. It is a skeleton of the catastrophe or skeleton of the mass-energy surface. The extremal point of this surface 
corresponds to the soliton solution. There exist also other types of solutions, with a different dynamical behavior. 
Perturbation theory gives small oscillations near the soliton solution. We went beyond perturbation theory, which 
helped us to investigate the stability of the star. 

Using the catastrophe theory, we can construct the mass-energy surface. Each section of this mass energy surface 
contains the different degrees of freedom of the star, which can be identified, by perturbation theory, with the 
characteristic frequencies. Using such sections, one can predict dynamical regimes of the star which cannot be 
described within the framework of perturbation theory. 

In order to predict the different regimes of behavior of the star, we construct the section M{l/R) (where R is the 
effective radius of the star) of the mass-energy surface M{R, N) at fixed A''. In other words, we construct an adiabatic 
potential [23]. The shape of M{l/R) follows from the bifurcation diagram. 

The type of regimes depends on the critical values , , , of the particle number. For the coupling constants 
a = {a/KW?) = 10 and /? = in the self-interacting potential (2.12), we find that Nc^ = 0.54, Nc^ = 0.68, A^Ca = 1-20 
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are the cuspoidal points, see Fig. 2. 



6.1. Stable soliton and oscillation 

From the bifurcation diagram (Fig. 2) it can be inferred that for N < Nc\ the function M{l/R) has only one 
minimum, which corresponds to the stable soliton solution (lower branch of the first cusp). 

The dependence of M on i? at some fixed value A'' < is schematically presented in Fig. 4. The smooth extremal 
point which corresponds to the minimum of the curve M{l/R) is associated with the soliton solution. The value of 
M in this point we call Mgouton- The marginal extremal point at i? = oo corresponds to the "homogeneous state", or 
plane wave solution in flat spacetime. In the self-generated gravitational field, it can be defined as an effectively free 
boson field solution for which the binding energy B — M — mN is vanishing. At infinite value of the radius R of the 
star the values M and N are final ones. This follows from the fact that in the limit R ^ oo the density of the star 
goes to zero (see Ref. [14], p. 29). The exact construction of the corresponding solution will be deferred to a future 
publication. These two extremal solution characterize a static configuration of the star. As we see from this figure, 
the "free boson" field is unstable. This extremum corresponds to the maximum. It means that the homogeneous state 
of the star will collapse from the size R = Rq = oo to the size R = Ri (see Fig. 4). From Ri, the size of the star 
will again increase. The value of R increases up to the initial value Rq. Thus the star will be in an oscillating regime. 
For such oscillations of the star, the curve presented in Fig. 4 is a kind of "adiabatic" potential which can also be 
obtained from a scale transformation. There are also other oscillations of the star in this potential that correspond 
to other values of M. For example, the other oscillation regime of the star corresponds to the horizontal line 1. In 
this case the amplitude of the oscillation is lower than the oscillations in the regime associated with the unstable free 
boson field described above. On the other hand, the oscillating regime, corresponding to the line 2, has a smaller 
amplitude than the one corresponding to the line 1 and so on. The limiting case is a stable soliton without oscillation. 
In the region of N < Nci, Fig. 4 gives a complete picture of the star's behavior. 

The interesting point here is that any arbitrary configuration of fixed particle number N cannot have a mass smaller 
than Msoiiton- The reason is that the configurations of the star are limited by the mapping F. For instance, the 
mapping F : {k,uj) — > {M,N) allows only a class of configurations away from the static soliton which preserve the 
integrals of motion M and N. Of course, the perturbed configurations of the stable soliton will evolve along the line 
dictated by the differential map. But there may arise other ways of choosing the parameters k and w. Moreover, one 
could think of extending the domain (fc, w) of the mapping F to a higher-dimensional space (fci, . . . , k^, Wi, . . . ,ujn), 
where fci, . . . , fcjv, wi, . . . , ojn correspond to additional degrees of freedom of the star. Such a space provides us with 
further lines of evolution of the perturbations, but the stability will still be determine by the bifurcation diagram 
M{N). Such a diagram comes from a numerical solution describing the stationary points of a Lagrange manifold of 
Einstein's equation. In this way one can describe virtually all perturbations which preserve the integrals of the motion 
M and A''. If we considered configurations which cannot be categorized by this choosen map, for example, when the 
total mass M and the total number of particles are not necessarily fixed, then one could presume that the star, 
via some oscillation process, settles down to a stationary stable configuration with some reduced mass Mreduced- 

The oscillations arise for M > Msoiuon of the mass of the star. In this region N < A^Cu at some value of 
M > Msoiiton, there may occur a gravitational collapse. We expect that the configuration will oscillate rather than 
collapse to a black hole not up to arbitrarily large mass M. There should exist the critical value of Mc = MschwarzschUd 
which depends on the radius of the star Rc with given mass Mc{Rc)- At this radius Rc the star will stop oscillating 
and start collapsing to a black hole. However, in order to show this rigorously, the analysis of Christodoulou [19] has 
to be extended to the case of the massive or even nonlinear scalar field. 

6.2. Collapse 

On the other hand, for A^ > Nc^, the section M{l/R) following from the bifurcation diagram of the mass-energy 
surface, has only the marginal extremum corresponding to the effectively free boson field solution. As we can infer from 
Fig. 5; the section M{l/R) is a montonically decreasing fimction. The homogeneous free field solution corresponds 
to the marginal extremum. It is a maximum of M{l/R). For the increase of the kinetic energy of the star, the radius 
decreases. This means that in "this region" of A'' there exists a collapse for any value of the mass. In the first stage, 
when Rstar » {ii/4:7T)Mschwarzschiid, it is a wave collapse (see [24] for details and references therein) which later 
induces the gravitational collapse. 

Due to this collapse, this state is unstable. The different horizontal lines, drawn in Fig. 5, correspond to different 
collapse regimes of the star with different initial radii. A similiar case is known for the creation of the two-dimensional 
plasma cavitons [25]. 
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6.3. Pulsation, oscillation, and collapse 

In the region < N < Nc^ the dependence of M on i? has as many extremal points as there exist branches 
of cusps at a given value of N. In the case where there are two cusp branches at given N, the dependence M{l/R) 
has two extremal points and a marginal extremum corresponding to the "free boson field" solution M = mN (Fig. 
6). The point of the minimum of M{l/R) corresponds to the stable soliton with the mass MgoUton- The instability 
of the homogeneous state results in an oscillation regime similar to that which has been described in Sect. 6.1. The 
maximum of M(l/R) corresponds to the unstable soliton with mass Mu- The instability of this soliton can occur in 
two possible ways. One of them is collapse, decreasing the radius of the star. The second one is increasing the star's 
radius (dispersion). Such an instability gives rise to an oscillation in the same manner as the marginal free particle 
extremum. There can be a lot of such oscillations (see Sect. 6.1) before the boson star will collapse. That means 
that this oscillating regime may be also unstable and, after several periods, collapses to the state of unstable soliton. 
Thus the existence of this unstable oscillating regime is due to the existence of the unstable soliton [24, 25]. 

Thus, in the most interesting region Nq^ < N < Nc^, the hmction M{l/R) has two extremal points, which 
correspond to the stable soliton (lower branch of the first cusp, sec Fig. 2) and the unstable soliton (higher branch 
of the first cusp), respectively, and one marginal extremum (free boson field solution). The shape of such a function 
indicates that at each value of N there exist the following configurations: a stable soliton with a mass Msoiiton{N) 
and unstable one with a mass Mu{N). At M{N) > Mu{N) it is very difficult to predict the evolution of the star. 

In the region Nci < N < Nc^ the dependence of M{l/R) may have more than two smooth extrema which are 
minima and maxima. For definiteness, let us consider the case when M{l/R) has three extremal points. Two of these 
points correspond to stable solitons with the mass Mg-^ and Ais2 ^-nd one corresponds to the unstable soliton with 
mass M„. There are two types of oscillations, corresponding to the minima Mg^ and Mg^ which exist at M < M„. 
In the first case, the maximal kinetic energy of the star corresponds to the radius of the first stable soliton. It is an 
oscillation in the first minimum. For the second case, the maximal kinetic energy corresponds to the radius of the 
second stable soliton. It is an oscillation in the second minimum. 

There is also a very interesting regime, for which the star has a mass M = Mu- In this case, there will exist a 
spatial type of oscillation in which it is difficult to predict in which direction (first or second minimum) the star will 
move from the point of maximum. It is connected with the indefiniteness of the behavior of the star in the state of 
the unstable soliton (dispersion or collapse). Such a regime of the star's behavior is characterized by a pulsation. The 
pulsation consists of at least two different types of oscillations. 

In the same manner one can analyze the case when M{l/R) has more than three extremal points. 

One general conclusion which can be drawn here is that in the region < N < there is a significant 

distribution of the instabilities of the star. The number of instabilities depends on the number of cusps and on the 
mutual branching of these cusps. If the next branch of some cusp corresponds to the higher values of the star's mass, 
then the number of instabilities increases by one. Vice versa, if the next branch of this cusp corresponds to lower 
values of the star's mass, then the number of instabilities decreases by one. 

7. Discussion 

Our stability criteria include the results of Refs. [10, 11, 12, 13] obtained by perturbation analysis. 

Until now, there exist some more or less successful attempts to prove the stability of the boson star. These results 
fight with the difficulty of the mathematical problem; they tried to solve the problem of stability in a quantitative 
way like Harrison et al. [8] or Shapiro and Teukolsky [7,9]. Gleiser [13] and Jetzer [12] got a upper limit for stability 
for the linear case and an additional | $ [^ potential. This limit was much higher than the first maximum in the 
(M,o-(0))-diagram. 

Later, Gleiser and Watkins [11] showed that in the linear case a change in stability occurs at the first extremum. At 

the following extrema, the higher modes (in the context of [8, 9]) are negative so that the star becomes more unstable. 
This result perfectly complies with those which are known from the analysis of neutron stars. Furthermore, Jetzer 
[12] showed that the zero node solution of the boson star is stable until one reaches the first cusp. Lee and Pang [10] 
did not require the particle number N to be constant, so they found that these solutions are unstable. 

The picture of the star's behavior, obtained here on the basis of the application of catastrophe theory to solitons 
[14], has a general character and has an analogous form for neutron stars [7; 8, see Fig. 8], fermion Q balls [26], and 
for dilaton stars [27]. In fact, for neutron stars the diagram M{N), obtained from numerical integration of a certain 
equation of state for cold, catalyzed matter, exhibit similar bifurcations as in the case of boson stars (sec Fig. 8 of 
Rcf. [8] and Fig. 47 of Ref. [7]). For the first cusp, Harrison ct al. [8] could even deduce correct stability criteria from 
the analysis of the mass-energy surface M{p,N), see Fig. 9 of Ref. [8]. Although this was done without knowledge 
of catastrophe theory or Arnold's classification of singularities, these 1965 results are in complete agreement with the 
more general criteria developed here. 

In comparison with Ref. [19] one should point out that in the hmit to ^ all solutions, corresponding to the stable 
solitons with < Mgoiiton < mN disappear and only unstable solutions remain. Such solutions disperse either to 
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infinity or collapse. The collapse configuration can form a black hole. This agrees completely with the results of Ref. 
[19]. 

After the submission of our paper, the recent paper of Scidcl and Sucn [28] appeared, in which the numerical evolu- 
tion of various configurations of boson stars using the full nonlinear Einstein equations have been studied. Moreover, 
perturbations which include a redistribution of scalar particles in the star and also accreation and annihilation of the 
bosons arc considered so that the total mass M and the total number of particles N arc not necessarily fixed. Their 
result, that the [/-branch star (unstable soliton) will either collapse to form a black hole or will disperse, agrees with 
our conclusion. The exception is that the unstable soliton will eventually settle down to a stable soliton. This is 
due to the possibility that the star is allowed to change the value of N or M . Such damping mechanisms have not 
been considered in our paper. The crucial role of the migration of unstable soliton ([/-branch) to a stable soliton 
(5-branch) has also been pointed out in Ref. [28]. According to this paper, a stable soliton (5-branch star), which 
is slightly perturbed, will oscillate with a fundamental frequency. This coincides with our conclusion about the oscil- 
lation regime of star near a stable soliton configuration. Thus the picture of the star's behavior obtained on the basis 
of theory of singularities of smooth maps (non -elementary catastrophe theory) completely coincides with conclusions 
obtained on the basis of extensive studies of the numerical evolution of boson stars. 
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FIG. 1. The Tolman mass M ( — ) in units of (l/mG) and particle number N ( ) in dimensionless units of (l/rn^G) as a 

function of the central density cr(0) = k/2 \ $(0) | for various a := (2a//tm^) and /3 = in the potential U. For a linear scalar 
field (a = 0) the Kaup limit Mxaup = 0.633 is recovered [17]. The maxima and minima correspond to the A2 singularities. 
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FIG. 2. The binding energy M — mN as a function of N at different values a — —5, 0, 5, 10 [17]. In this bifurcation diagram, 
the lower branch of each cusp corresponds to a stable star configuration. See, for comparison. Fig. 1. 
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FIG. 3. "Magnified" view of the binding energy for the same parameters as in Fig. 
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FIG. 4. "Adiabatic" potential M = M(l/R) of the star for < Nc^ (Schematic construction following Ref. [22, 24]). 
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FIG. 5. "Adiabatic" potential M = M(l/R) of the star for > Nc^, describing the collapse. 
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FIG. 6. "Adiabatic" potential M 



= M{l/R) of the star for A'' e [Nc2,Nc3], describing oscillation and collapse. 
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